Übung Multivariate Verfahren
Janika Saretzki
2025-01-23
Allgemeine Anmerkungen
#: Wird im R-Markdown-Text für Überschriften verwendet. Mehrere Ebenen können durch Hinzufügen weiterer Rauten erstellt werden.**Wort**: Wort wird dick gedruckt dargestellt.*Wort*: Wort wird kursiv dargestellt.Sonderzeichen können durch einen Backslash (
\) “escapet” werden, damit sie als normale Zeichen erscheinen.
Lineare Gemischte Modelle
Beispiel 1
Datensatz “LMM.xlsx”
# install.packages("readxl")
# library(readxl)
lmm = read_excel("LMM.xlsx")
head(lmm) # Zeigt erste Zeilen des Datensatzes "lmm"## # A tibble: 6 × 4
## ids department experience salary
## <dbl> <chr> <dbl> <dbl>
## 1 1 sociology 0 36095.
## 2 2 biology 1 55254.
## 3 3 english 1 59140.
## 4 4 informatics 8 78325.
## 5 5 statistics 1 74054.
## 6 6 sociology 3 46508.summary(lmm) # Statistische Zusammenfassung des Datensatzes "lmm" einschließlich Minimum, Median, Mittelwert, Maximum und Quartile für numerische Variablen## ids department experience salary
## Min. : 1.00 Length:100 Min. :0.00 Min. :36095
## 1st Qu.: 25.75 Class :character 1st Qu.:2.00 1st Qu.:52125
## Median : 50.50 Mode :character Median :4.00 Median :62229
## Mean : 50.50 Mean :4.39 Mean :64278
## 3rd Qu.: 75.25 3rd Qu.:7.00 3rd Qu.:78187
## Max. :100.00 Max. :9.00 Max. :90027# View(lmm) # Öffnet tabellarische Ansicht des Datensatzes "lmm" im RStudio Viewer
str(lmm) # Zeigt Struktur des Datensatzes "lmm" einschließlich Variablentypen und ersten Beispieldaten)## tibble [100 × 4] (S3: tbl_df/tbl/data.frame)
## $ ids : num [1:100] 1 2 3 4 5 6 7 8 9 10 ...
## $ department: chr [1:100] "sociology" "biology" "english" "informatics" ...
## $ experience: num [1:100] 0 1 1 8 1 3 8 2 3 2 ...
## $ salary : num [1:100] 36095 55254 59140 78325 74054 ...names(lmm) # Gibt Namen der Variablen (Spalten) des Datensatzes "lmm" aus## [1] "ids" "department" "experience" "salary"unique(lmm$department) # Zeigt alle einzigartigen Werte der Variable "department"## [1] "sociology" "biology" "english" "informatics" "statistics"lapply(lmm[, c("department")], unique) # Funktion "lappy()": Wendet ausgewählte Funktion (z.B. "unique()") auf mehrere Variablen des Datensatzes an.## $department
## [1] "sociology" "biology" "english" "informatics" "statistics"- Hierarchische Struktur der Daten: Jeder Dozent bzw. Dozentin (id) ist einer Abteilung (department) zugeordnet.
- Annahme: Die Abteilungen haben einen wesentlichen Einfluss auf das Gehalt (salary).
Einfaches Lineares Modell
\(\rightarrow\) Einfaches Modell zur Schätzung durchschnittlicher Gehaltsunterschiede zwischen Abteilungen mittels fester Effekte (Fixed Effects Model)
model1 = lm(salary ~ department, lmm)
summary(model1)##
## Call:
## lm(formula = salary ~ department, data = lmm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11922.3 -2985.0 -52.7 3150.7 11138.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51843 1127 45.983 < 0.0000000000000002 ***
## departmentenglish 10617 1594 6.659 0.00000000178 ***
## departmentinformatics 26080 1594 16.357 < 0.0000000000000002 ***
## departmentsociology -3826 1594 -2.399 0.0184 *
## departmentstatistics 29304 1594 18.379 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5042 on 95 degrees of freedom
## Multiple R-squared: 0.8809, Adjusted R-squared: 0.8759
## F-statistic: 175.6 on 4 and 95 DF, p-value: < 0.00000000000000022qqnorm(resid(model1))
shapiro.test(resid(model1))##
## Shapiro-Wilk normality test
##
## data: resid(model1)
## W = 0.99367, p-value = 0.925hist(resid(model1))
Unconditional Random Effects Model (UREM)
- Flexiblere Modellierung der Gehaltsunterschiede zwischen den Abteilungen als zufällige Effekte
- Random Intercept Modell ohne Prädiktor)
- Annahme: Jede Abteilung hat ein eigenes durchschnittliches Gehaltsniveau, das um den globalen Mittelwert variiert
- Nutzen: Isolierung der reinen Gruppeneffekte als Grundlage für die spätere Einbeziehung individueller Prädiktoren [z.B. Berufserfahrung (experience)]
# install.packages("lme4")
# library(lme4)
# install.packages("lmerTest")
# library(lmerTest)
model2 = lmer(salary ~ 1 + (1|department), lmm)
summary(model2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ 1 + (1 | department)
## Data: lmm
##
## REML criterion at convergence: 1994.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.38295 -0.58402 -0.00795 0.61588 2.22449
##
## Random effects:
## Groups Name Variance Std.Dev.
## department (Intercept) 221997992 14900
## Residual 25422058 5042
## Number of obs: 100, groups: department, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 64278 6682 4 9.619 0.000653 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(model2)$coefficients # $coefficients um nur Fixed Effects zu zeigen## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 64277.74 6682.351 4 9.61903 0.0006530776VarCorr(model2) # Rohstruktur der Varianz-Kovarianz Matrix der Random Effects, beinhaltet lediglich Standardabweichungen ## Groups Name Std.Dev.
## department (Intercept) 14900
## Residual 5042as.data.frame(VarCorr(model2)) # Legt vollständige Berechnungen (Varianz UND Standardabweichung) offen## grp var1 var2 vcov sdcor
## 1 department (Intercept) <NA> 221997992 14899.597
## 2 Residual <NA> <NA> 25422058 5042.029confint(model2, level = 0.95) # Output-Erweiterung: Konfidenzintervalle## Computing profile confidence intervals ...## 2.5 % 97.5 %
## .sig01 7911.713 29098.712
## .sigma 4401.953 5853.681
## (Intercept) 49907.902 78647.572Intra-Klassen-Korrelation (engl. Intra-Class-Correlation; ICC)
- ICC berechnet den Anteil der Gesamtvarianz, der durch Unterschiede zwischen den Abteilungen erklärt wird
- ICC = Varianz zwischen den Gruppen / Gesamtvarianz
ICC_model2 = 221997992 / (221997992 + 25422058)
ICC_model2## [1] 0.8972514- 90% der Gesamtvarianz der Gehälter wird durch Unterschiede zwischen den Abteilungen erklärt
- Restliche 10%? Erklärt durch andere Variablen (z.B. Unterschiede auf Individualebene, d.h. zwischen den Personen) und Fehler
Random Intercept Fixed Slope Model
- Erweiterung der Analyse, um individuelle Merkmale [z.B. Berufserfahrung (experience)] als Prädiktoren zu berücksichtigen
- Annahme: Berufserfahrung (experience) hat einen linearen Einfluss auf das Gehalt, der über alle Abteilungen hinweg konstant bleibt (Fixed Slope)
- Nutzen: Kombination von individuellen Effekten (experience) und Gruppeneffekten (Abteilungsunterschiede, Random Intercept) zur differenzierten Analyse der Gehaltsunterschiede
model3 = lmer(salary ~ experience + (1|department), lmm)
summary(model3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ experience + (1 | department)
## Data: lmm
##
## REML criterion at convergence: 1953.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.93140 -0.80670 -0.07917 0.80156 2.12467
##
## Random effects:
## Groups Name Variance Std.Dev.
## department (Intercept) 215356636 14675
## Residual 18880965 4345
## Number of obs: 100, groups: department, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 60469.557 6609.551 4.079 9.149 0.000723 ***
## experience 867.468 148.681 94.013 5.834 0.0000000762 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## experience -0.099# REML durch ML ersetzen ...?
# model3 = lmer(salary ~ experience + (1|department), REML = FALSE, lmm)
# summary(model3)
summary(model3)$coefficients ## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 60469.5567 6609.5511 4.07891 9.148814 0.00072262119073
## experience 867.4677 148.6809 94.01310 5.834426 0.00000007616036confint(model3, level = 0.95)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## .sig01 7805.4756 28645.859
## .sigma 3773.5921 5018.089
## (Intercept) 46288.9712 74656.162
## experience 575.0961 1160.768sqrt(221997992)## [1] 14899.6# install.packages("ggplot2")
# library(ggplot2)
lmm$random.intercept.fixed.slope = predict(model3)
ggplot(lmm, aes(x = experience, y = salary, color = department)) +
geom_point() +
geom_line(aes(y = random.intercept.fixed.slope)) +
labs(x = "Experience", y = "Salary", color = "Department") +
theme_minimal() +
theme(panel.grid = element_blank(), axis.line = element_line())
Random Intercept Random Slope Model
- Erweiterung: Modellierung sowohl der mittleren Gehälter (Random Intercepts) als auch des Einflusses der Berufserfahrung (Random Slopes) als abteilungsabhängige Effekte
- Annahme: Die mittleren Gehälter variieren zwischen den Abteilungen UND der Einfluss der Berufserfahrung auf das Gehalt unterscheidet sich ebenfalls zwischen den Abteilungen
model4 = lmer(salary ~ experience + (experience|department), lmm)
summary(model4)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ experience + (experience | department)
## Data: lmm
##
## REML criterion at convergence: 1941.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76050 -0.84375 -0.01161 0.73272 2.32260
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## department (Intercept) 227233147 15074.3
## experience 469881 685.5 -0.26
## Residual 15518494 3939.4
## Number of obs: 100, groups: department, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 60149.368 6779.433 4.000 8.872 0.000892 ***
## experience 922.801 335.170 3.987 2.753 0.051390 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## experience -0.274summary(model4)$coefficients ## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 60149.3685 6779.4330 3.999646 8.872330 0.0008917901
## experience 922.8008 335.1701 3.986902 2.753231 0.0513899885confint(model4, level = 0.95)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## .sig01 7944.7287779 29496.9246374
## .sig02 -0.8478171 0.6389808
## .sig03 256.5642789 1430.1370444
## .sigma 3427.0528130 4592.9136408
## (Intercept) 45594.6492882 74748.1610306
## experience 206.2522390 1648.7374105lmm$random.intercept.random.slope = predict(model4)
ggplot(lmm, aes(x = experience, y = salary, color = department)) +
geom_point() +
geom_line(aes(y = random.intercept.random.slope)) +
labs(x = "Experience", y = "Salary", color = "Department") +
theme_minimal() +
theme(panel.grid = element_blank(), axis.line = element_line())
Chi-Quadrat-Likelihood Ratio Test
# Model 3: Random Intercept, Fixed Slope
# Model 4: Random Intercept, Random Slope
anova(model3, model4)## refitting model(s) with ML (instead of REML)## Data: lmm
## Models:
## model3: salary ~ experience + (1 | department)
## model4: salary ~ experience + (experience | department)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## model3 4 1992.2 2002.6 -992.10 1984.2
## model4 6 1986.5 2002.1 -987.23 1974.5 9.7504 2 0.007634 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Fixed Intercept Random Slope Model
- Der Intercept wird für alle Abteilungen als gleich angenommen (Fixed Intercept)
- Der Einfluss der Berufserfahrung darf zwischen den Abteilungen variieren (Random Slope)
Warum ist dieses Modell nicht sinnvoll?
- Aus dem Random Intercept Fixed Slope Modell wissen wir, dass sich die Abteilungen im durchschnittlichen Gehalt (Intercept) unterscheiden. Die Annahme eines festen Intercepts für alle Abteilungen ist daher nicht realistisch.
- Wenn die mittleren Gehälter der Abteilungen unterschiedlich sind, würde ein Modell mit festem Intercept wichtige Variabilität in den Daten ignorieren. Das Modell vernachlässigt daher die Heterogenität der Abteilungen hinsichtlich der mittleren Gehälter.
model5 = lmer(salary ~ experience + (0 + experience | department), lmm)
summary(model5)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: salary ~ experience + (0 + experience | department)
## Data: lmm
##
## REML criterion at convergence: 2080.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.67330 -0.57808 -0.03245 0.60669 2.86364
##
## Random effects:
## Groups Name Variance Std.Dev.
## department experience 5309875 2304
## Residual 74068552 8606
## Number of obs: 100, groups: department, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 59102.075 1546.912 94.013 38.206 <0.0000000000000002 ***
## experience 1065.242 1071.503 4.445 0.994 0.371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## experience -0.227lmm$fixed.intercept.random.slope = predict(model5)
ggplot(lmm, aes(x = experience, y = salary, color = department)) +
geom_point() +
geom_line(aes(y = fixed.intercept.random.slope)) +
labs(x = "Experience", y = "Salary", color = "Department") +
theme_minimal() +
theme(panel.grid = element_blank(), axis.line = element_line())
Beispiel 2
Datensatz “LMM_2.xlsx”
# install.packages("readxl")
# library(readxl)
performance_initial = read_excel("LMM_2.xlsx")
names(performance_initial)## [1] "subject" "trainer" "training"
## [4] "performance_time0" "performance_time1" "performance_time2"
## [7] "performance_time3" "performance_time4" "performance_time5"
## [10] "performance_time6" "performance_time7" "performance_time8"
## [13] "performance_time9" "performance_time10"str(performance_initial)## tibble [80 × 14] (S3: tbl_df/tbl/data.frame)
## $ subject : num [1:80] 1 2 3 4 5 6 7 8 9 10 ...
## $ trainer : num [1:80] 5 5 5 5 5 5 5 5 5 5 ...
## $ training : chr [1:80] "control" "control" "control" "control" ...
## $ performance_time0 : num [1:80] 38.9 32.9 35.3 39.3 36.8 ...
## $ performance_time1 : num [1:80] 37 36.4 36.9 41.7 34.3 ...
## $ performance_time2 : num [1:80] 37.9 31.7 37.1 41.4 35.5 ...
## $ performance_time3 : num [1:80] 37.6 34 36.4 41.1 33 ...
## $ performance_time4 : num [1:80] 38 33.4 36.5 40.3 36.6 ...
## $ performance_time5 : num [1:80] 37.2 32.6 34.9 41.1 38.3 ...
## $ performance_time6 : num [1:80] 37.1 33.2 38 40.4 35.3 ...
## $ performance_time7 : num [1:80] 37.3 35.5 39.6 39.1 33.5 ...
## $ performance_time8 : num [1:80] 34.5 34.8 40 36.2 33.7 ...
## $ performance_time9 : num [1:80] 37.8 32.3 37.3 39.3 35.3 ...
## $ performance_time10: num [1:80] 34.5 32.8 40.9 38.2 35.9 ...lapply(performance_initial[, c("training")], unique) # EXTRAKTION EINDEUTIGER WERTE EINER AUSWAHL AN VARIABLEN (HIER: VARIABLE "training" AUS DATENSATZ "LMM_2")## $training
## [1] "control" "training"head(performance_initial)## # A tibble: 6 × 14
## subject trainer training performance_time0 performance_time1 performance_time2
## <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 1 5 control 38.9 37.0 37.9
## 2 2 5 control 32.9 36.4 31.7
## 3 3 5 control 35.3 36.9 37.1
## 4 4 5 control 39.3 41.7 41.4
## 5 5 5 control 36.8 34.3 35.5
## 6 6 5 control 43.6 42.3 43.5
## # ℹ 8 more variables: performance_time3 <dbl>, performance_time4 <dbl>,
## # performance_time5 <dbl>, performance_time6 <dbl>, performance_time7 <dbl>,
## # performance_time8 <dbl>, performance_time9 <dbl>, performance_time10 <dbl># WIDE-FORMAT >> LONG-FORMAT
# install.packages("tidyverse") # ÜBERTRAGEN IN SETUP-CHUNK
# library(tidyverse) #ÜBERTRAGEN IN SETUP-CHUNK
performance = pivot_longer(
data = performance_initial,
cols = c(performance_time0, performance_time1, performance_time2,
performance_time3, performance_time4, performance_time5,
performance_time6, performance_time7, performance_time8,
performance_time9, performance_time10),
names_to = "Time", # Der Name der ursprünglichen Spalten (z. B. "performance_time0") wird in eine neue Variable namens "Time" geschrieben.
names_prefix = "performance_time", # Der Präfix "performance_time" wird aus den Spaltennamen entfernt, sodass nur die Zeit (z. B. "0", "1", ...) übrig bleibt.
values_to = "Performance" # Die Werte in den ausgewählten Spalten werden in eine neue Spalte namens "Performance" geschrieben.
)
performance$Time = as.numeric(performance$Time)
names(performance)## [1] "subject" "trainer" "training" "Time" "Performance"head(performance)## # A tibble: 6 × 5
## subject trainer training Time Performance
## <dbl> <dbl> <chr> <dbl> <dbl>
## 1 1 5 control 0 38.9
## 2 1 5 control 1 37.0
## 3 1 5 control 2 37.9
## 4 1 5 control 3 37.6
## 5 1 5 control 4 38.0
## 6 1 5 control 5 37.2# install.packages("writexl")
# library(writexl)
writexl::write_xlsx(performance, "LMM_2_Long_Format.xlsx")Die Daten haben eine Längsschnittstruktur, da für jede Person (subject) die Outcome-Variable (Performance, d.h. Sprintzeit) zu mehreren Zeitpunkten (Time) erfasst wurde.
\(\rightarrow\) Diese Struktur ermöglicht es uns, die individuelle Entwicklung der Sprintleistung im Zeitverlauf zu analysieren.
Zusätzlich enthalten die Daten zwei Gruppierungsvariablen:
Trainer
Trainingsstatus
Mögliche Fragestellungen: Wie verändert sich die Sprintzeit im Durchschnitt über die Zeit in Abhängigkeit von der Teilnahme an einem intensiven Intervalltraining?
Mögliche Hypothese: Athleten in der Trainingsgruppe zeigen im Durchschnitt eine stärkere Abnahme der Sprintzeit über die Zeit als Athleten in der Kontrollgruppe.
Einfaches Wachstumskurvenmodell (ohne Prädiktor)
Random Intercept Fixed Slope Model
Annahme: Jeder Athlet eine individuelle Ausgangsleistung (Intercept, d.h. anfängliche Sprintzeit) hat, die um einen globalen Mittelwert streut, während die Wachstumsrate (Steigung, d.h. Verbesserung der Sprintzeit) für alle Athleten gleich bleibt.
\(\rightarrow\) Einfacher Einstieg, da es individuelle Unterschiede im Ausgangspunkt (Intercept) berücksichtigt, ohne die Komplexität zusätzlicher Prädiktoren wie Trainingsstatus oder Trainer oder deren Interaktion einzubeziehen. Es erlaubt eine erste Einschätzung der zeitlichen Entwicklung der Sprintzeiten und der Variabilität zwischen den Athleten.
model6 = lmer(Performance ~ Time + (1|subject), performance)
summary(model6)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Performance ~ Time + (1 | subject)
## Data: performance
##
## REML criterion at convergence: 3673.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1338 -0.6552 0.0057 0.6154 3.1002
##
## Random effects:
## Groups Name Variance Std.Dev.
## subject (Intercept) 6.683 2.585
## Residual 2.804 1.675
## Number of obs: 880, groups: subject, 80
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 37.47091 0.30771 94.10418 121.77 <0.0000000000000002 ***
## Time -0.27717 0.01785 799.00000 -15.53 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Time -0.290performance$Random.Intercept.Fixed.Slope = predict(model6, re.form = NA)
# Zusatzargument "re.form = NA" - GLOBALE STEIGUNG, keine individuelle STEIGUNG
ggplot(performance, aes(x = Time, y = Performance, group = subject)) +
geom_point(size = 0.1) +
geom_line(size = 0.1) +
geom_line(aes(x = Time, y = Random.Intercept.Fixed.Slope), color = "black", size = 0.8) +
labs(x = "Time", y = "Performance") +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.line = element_line(),
legend.position = "none"
)
Random Intercept Random Slope Model
model7 = lmer(Performance ~ Time + (Time|subject), performance)
summary(model7)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Performance ~ Time + (Time | subject)
## Data: performance
##
## REML criterion at convergence: 3559.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.15901 -0.64330 -0.02121 0.59911 2.98766
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 9.13058 3.0217
## Time 0.06458 0.2541 -0.52
## Residual 2.10207 1.4499
## Number of obs: 880, groups: subject, 80
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 37.47091 0.34999 78.99982 107.06 < 0.0000000000000002 ***
## Time -0.27717 0.03234 78.99935 -8.57 0.000000000000674 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Time -0.547Konditionales Wachstumskurvenmodell (mit Prädiktor)
Random Intercept Random Slope Model (2 Gruppen)
- Zusätzlicher Prädiktor: Beeinflusst die Teilnahme an einem Training die Sprintzeit über die Zeit?
Das Modell berücksichtigt:
Random Intercept: Jede Person (subject) hat einen individuellen Ausgangswert (Intercept), der um den globalen Mittelwert streut.
Random Slope: Jede Person hat eine individuelle Wachstumsrate (Slope), die um die globale Steigung streut.
Training als fester Prädiktor: Das Modell testet, ob die Sprintzeit im Zeitverlauf durch die Teilnahme am Training beeinflusst wird. Die Variable Training hat zwei Ausprägungen: Training vs. Control
model8 = lmer(Performance ~ Time * training + (Time|subject), performance)
summary(model8)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Performance ~ Time * training + (Time | subject)
## Data: performance
##
## REML criterion at convergence: 3555.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.13461 -0.65993 0.00125 0.61249 2.96479
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subject (Intercept) 9.21415 3.0355
## Time 0.05937 0.2437 -0.53
## Residual 2.10208 1.4499
## Number of obs: 880, groups: subject, 80
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 37.26848 0.49707 78.00008 74.977 < 0.0000000000000002
## Time -0.19897 0.04430 78.00246 -4.492 0.0000241
## trainingtraining 0.40486 0.70296 78.00008 0.576 0.5663
## Time:trainingtraining -0.15639 0.06264 78.00246 -2.497 0.0146
##
## (Intercept) ***
## Time ***
## trainingtraining
## Time:trainingtraining *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Time trnngt
## Time -0.551
## tranngtrnng -0.707 0.390
## Tm:trnngtrn 0.390 -0.707 -0.551performance$Random.Intercept.Random.Slope = predict(model8, re.form = NA)
ggplot(performance, aes(x = Time, y = Performance, group = subject)) +
geom_point(aes(color = training), size = 0.1) +
geom_line(aes(color = training), size = 0.1) +
geom_line(aes(x = Time, y = Random.Intercept.Random.Slope), color = "black", size = 1) +
stat_summary(fun = mean, aes(x = Time, y = Performance, group = training, color = training), geom = "line", linetype = "dashed", size = 1) +
labs(x = "Time", y = "Performance", color = "Group") +
theme_minimal() +
theme(panel.grid = element_blank(), axis.line = element_line(), legend.position = "bottom"
)
Explorative Faktorenanalyse (EFA)
Fragestellung: Auf wie viele Faktoren lassen sich die Zusammenhänge zwischen den gegebenen neun Skalen zurückführen?
Stichprobenbeschreibung
nrow(EFA) # Anzahl der Fälle insgesamt## [1] 299# Geschlecht
table(EFA$SEX)##
## 1 2
## 139 160str(EFA$SEX)## dbl+lbl [1:299] 2, 2, 2, 1, 1, 1, 2, 2, 2, 2, 1, 2, 2, 1, 2, 2, 1, 2, 2, 2...
## @ label : chr "Geschlecht"
## @ format.spss : chr "F1.0"
## @ display_width: int 5
## @ labels : Named num [1:2] 1 2
## ..- attr(*, "names")= chr [1:2] "männlich" "weiblich"EFA$SEX = as.factor(EFA$SEX)
levels(EFA$SEX) = c("männlich", "weiblich")
table(EFA$SEX)##
## männlich weiblich
## 139 160prop.table(table(EFA$SEX)) * 100##
## männlich weiblich
## 46.48829 53.51171round(prop.table(table(EFA$SEX)),3) * 100##
## männlich weiblich
## 46.5 53.5# Alter
min(EFA$AGE)## [1] 15max(EFA$AGE)## [1] 81# EFA = EFA[EFA$AGE >= 18 & EFA$AGE <= 65, ]
EFA = subset(EFA, select = -c(SEX, AGE))
names(EFA)## [1] "Risikobereitschaft" "Impulsivität"
## [3] "Verantwortungslosigkeit" "Aktivität"
## [5] "Kontaktfreudigkeit" "Selbstbewusstsein"
## [7] "Minderwertigkeit" "Schwermut"
## [9] "Besorgtheit"EFA = EFA %>%
rename(
V1 = Risikobereitschaft,
V2 = Impulsivität,
V3 = Verantwortungslosigkeit,
V4 = Aktivität,
V5 = Kontaktfreudigkeit,
V6 = Selbstbewusstsein,
V7 = Minderwertigkeit,
V8 = Schwermut,
V9 = Besorgtheit
)
names(EFA)## [1] "V1" "V2" "V3" "V4" "V5" "V6" "V7" "V8" "V9"cor(EFA)## V1 V2 V3 V4 V5 V6
## V1 1.00000000 0.47351664 0.44323461 0.1711286 0.26041822 0.27754924
## V2 0.47351664 1.00000000 0.39993072 0.2024946 0.26048243 0.15869634
## V3 0.44323461 0.39993072 1.00000000 -0.3245046 0.07312967 0.07172468
## V4 0.17112859 0.20249461 -0.32450459 1.0000000 0.45186569 0.33243051
## V5 0.26041822 0.26048243 0.07312967 0.4518657 1.00000000 0.39587257
## V6 0.27754924 0.15869634 0.07172468 0.3324305 0.39587257 1.00000000
## V7 -0.19405554 0.06794166 0.08302704 -0.3622716 -0.40195207 -0.47485579
## V8 -0.07104343 0.13201378 0.19287894 -0.3474677 -0.31153103 -0.25263399
## V9 -0.17214201 0.20128015 0.11456815 -0.1665212 -0.24523910 -0.25264940
## V7 V8 V9
## V1 -0.19405554 -0.07104343 -0.1721420
## V2 0.06794166 0.13201378 0.2012802
## V3 0.08302704 0.19287894 0.1145681
## V4 -0.36227155 -0.34746771 -0.1665212
## V5 -0.40195207 -0.31153103 -0.2452391
## V6 -0.47485579 -0.25263399 -0.2526494
## V7 1.00000000 0.71410876 0.6765718
## V8 0.71410876 1.00000000 0.7026063
## V9 0.67657185 0.70260633 1.0000000cor(EFA[, c("V1", "V2", "V3", "V4", "V5", "V6")])## V1 V2 V3 V4 V5 V6
## V1 1.0000000 0.4735166 0.44323461 0.1711286 0.26041822 0.27754924
## V2 0.4735166 1.0000000 0.39993072 0.2024946 0.26048243 0.15869634
## V3 0.4432346 0.3999307 1.00000000 -0.3245046 0.07312967 0.07172468
## V4 0.1711286 0.2024946 -0.32450459 1.0000000 0.45186569 0.33243051
## V5 0.2604182 0.2604824 0.07312967 0.4518657 1.00000000 0.39587257
## V6 0.2775492 0.1586963 0.07172468 0.3324305 0.39587257 1.00000000cor(EFA[, grep("V[1-9]",names(EFA))])## V1 V2 V3 V4 V5 V6
## V1 1.00000000 0.47351664 0.44323461 0.1711286 0.26041822 0.27754924
## V2 0.47351664 1.00000000 0.39993072 0.2024946 0.26048243 0.15869634
## V3 0.44323461 0.39993072 1.00000000 -0.3245046 0.07312967 0.07172468
## V4 0.17112859 0.20249461 -0.32450459 1.0000000 0.45186569 0.33243051
## V5 0.26041822 0.26048243 0.07312967 0.4518657 1.00000000 0.39587257
## V6 0.27754924 0.15869634 0.07172468 0.3324305 0.39587257 1.00000000
## V7 -0.19405554 0.06794166 0.08302704 -0.3622716 -0.40195207 -0.47485579
## V8 -0.07104343 0.13201378 0.19287894 -0.3474677 -0.31153103 -0.25263399
## V9 -0.17214201 0.20128015 0.11456815 -0.1665212 -0.24523910 -0.25264940
## V7 V8 V9
## V1 -0.19405554 -0.07104343 -0.1721420
## V2 0.06794166 0.13201378 0.2012802
## V3 0.08302704 0.19287894 0.1145681
## V4 -0.36227155 -0.34746771 -0.1665212
## V5 -0.40195207 -0.31153103 -0.2452391
## V6 -0.47485579 -0.25263399 -0.2526494
## V7 1.00000000 0.71410876 0.6765718
## V8 0.71410876 1.00000000 0.7026063
## V9 0.67657185 0.70260633 1.0000000cor(EFA[, grep("^V[1-9]$",names(EFA))])## V1 V2 V3 V4 V5 V6
## V1 1.00000000 0.47351664 0.44323461 0.1711286 0.26041822 0.27754924
## V2 0.47351664 1.00000000 0.39993072 0.2024946 0.26048243 0.15869634
## V3 0.44323461 0.39993072 1.00000000 -0.3245046 0.07312967 0.07172468
## V4 0.17112859 0.20249461 -0.32450459 1.0000000 0.45186569 0.33243051
## V5 0.26041822 0.26048243 0.07312967 0.4518657 1.00000000 0.39587257
## V6 0.27754924 0.15869634 0.07172468 0.3324305 0.39587257 1.00000000
## V7 -0.19405554 0.06794166 0.08302704 -0.3622716 -0.40195207 -0.47485579
## V8 -0.07104343 0.13201378 0.19287894 -0.3474677 -0.31153103 -0.25263399
## V9 -0.17214201 0.20128015 0.11456815 -0.1665212 -0.24523910 -0.25264940
## V7 V8 V9
## V1 -0.19405554 -0.07104343 -0.1721420
## V2 0.06794166 0.13201378 0.2012802
## V3 0.08302704 0.19287894 0.1145681
## V4 -0.36227155 -0.34746771 -0.1665212
## V5 -0.40195207 -0.31153103 -0.2452391
## V6 -0.47485579 -0.25263399 -0.2526494
## V7 1.00000000 0.71410876 0.6765718
## V8 0.71410876 1.00000000 0.7026063
## V9 0.67657185 0.70260633 1.0000000round(cor(EFA), 2)## V1 V2 V3 V4 V5 V6 V7 V8 V9
## V1 1.00 0.47 0.44 0.17 0.26 0.28 -0.19 -0.07 -0.17
## V2 0.47 1.00 0.40 0.20 0.26 0.16 0.07 0.13 0.20
## V3 0.44 0.40 1.00 -0.32 0.07 0.07 0.08 0.19 0.11
## V4 0.17 0.20 -0.32 1.00 0.45 0.33 -0.36 -0.35 -0.17
## V5 0.26 0.26 0.07 0.45 1.00 0.40 -0.40 -0.31 -0.25
## V6 0.28 0.16 0.07 0.33 0.40 1.00 -0.47 -0.25 -0.25
## V7 -0.19 0.07 0.08 -0.36 -0.40 -0.47 1.00 0.71 0.68
## V8 -0.07 0.13 0.19 -0.35 -0.31 -0.25 0.71 1.00 0.70
## V9 -0.17 0.20 0.11 -0.17 -0.25 -0.25 0.68 0.70 1.00# install.packages("corrplot")
# library(corrplot)
corrplot(cor(EFA), method = "color")
# colorbrewer2.org
corrplot(cor(EFA), method = "color", col = colorRampPalette(c("#c51b7d", "white", "#4d9221"))(200))
Anwendungsvoraussetzungen
# install.packages("psych")
# library(psych)
# KMO Kriterium (Kaiser-Meyer-Olkin Kriterium)
KMO(EFA) # MSA = Measure of Sampling Adequacy## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = EFA)
## Overall MSA = 0.71
## MSA for each item =
## V1 V2 V3 V4 V5 V6 V7 V8 V9
## 0.65 0.66 0.51 0.59 0.83 0.80 0.79 0.76 0.73# Bartlett-Test auf Sphärizität
# H0: Korrelationsmatrix ist Einheitsmatrix
cortest.bartlett(cor(EFA))## $chisq
## [1] 343.2489
##
## $p.value
## [1] 0.0000000000000000000000000000000000000000000000000008835381
##
## $df
## [1] 36Bestimmung der Anzahl der Faktoren
# Zugang über EIGENWERTE
# KAISER KRITERIUM
## Nur Faktorenwerte mit einem Eigenwert > 1 werden extrahiert
eigen(cor(EFA))## eigen() decomposition
## $values
## [1] 3.1929551 2.0686210 1.1907609 0.7097815 0.5906573 0.4437106 0.3270208
## [8] 0.2445874 0.2319054
##
## $vectors
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.20180955 -0.48942333 0.20344713 0.26857384 -0.41978643 -0.57276092
## [2,] 0.05295321 -0.56382622 -0.19165245 0.30484118 -0.04053450 0.61663804
## [3,] -0.05434755 -0.50561927 0.48376815 -0.11450019 0.27195504 0.11761254
## [4,] 0.33517591 -0.02651339 -0.61858851 0.23814739 -0.20414333 -0.07730031
## [5,] 0.35670852 -0.20780049 -0.26760594 -0.09271811 0.77058157 -0.32603363
## [6,] 0.34529673 -0.17439342 -0.09410667 -0.80310917 -0.32503994 0.10879561
## [7,] -0.48333867 -0.11258204 -0.18216456 0.07951369 0.02199774 -0.18949111
## [8,] -0.44015227 -0.23338639 -0.18078765 -0.26969617 -0.07327041 -0.32769616
## [9,] -0.40908585 -0.21324391 -0.39857012 -0.18021110 0.03468271 0.10231488
## [,7] [,8] [,9]
## [1,] -0.04999693 -0.18684477 0.25402497
## [2,] -0.39873981 -0.04236226 -0.08344843
## [3,] 0.54950545 0.27218010 -0.17442216
## [4,] 0.49054791 0.28578525 -0.27861790
## [5,] -0.19628183 -0.05439266 0.08825754
## [6,] -0.12730749 0.20320571 0.14672936
## [7,] -0.19877833 0.72621308 0.33303696
## [8,] -0.18582194 -0.16396396 -0.68717485
## [9,] 0.40882395 -0.45718542 0.45608648eigen(cor(EFA))$values## [1] 3.1929551 2.0686210 1.1907609 0.7097815 0.5906573 0.4437106 0.3270208
## [8] 0.2445874 0.2319054# SCREE Plot
psych::scree(cor(EFA))
# PARALLELANALYSE
fa.parallel(EFA, fm="ml")
## Parallel analysis suggests that the number of factors = 3 and the number of components = 3Faktorenrotation
fa(EFA, nfactors = 3, fm = "ml", rotate = "promax")$loadings## Loading required namespace: GPArotation##
## Loadings:
## ML1 ML3 ML2
## V1 -0.101 0.672
## V2 0.331 0.617 0.266
## V3 0.864 -0.494
## V4 -0.230 1.032
## V5 -0.189 0.209 0.395
## V6 -0.280 0.217 0.233
## V7 0.825
## V8 0.817 0.162
## V9 0.913 0.176
##
## ML1 ML3 ML2
## SS loadings 2.426 1.752 1.632
## Proportion Var 0.270 0.195 0.181
## Cumulative Var 0.270 0.464 0.646# Promax: Oblique Rotationsmethode; Faktoren sind nicht unkorreliert, d.h. nicht unabhängig voneinanderKommunalitäten
fa(EFA, nfactors = 3, fm = "ml", rotate = "promax")$communalities## V1 V2 V3 V4 V5 V6 V7 V8
## 0.5042159 0.5318089 0.6699846 0.8576953 0.3831535 0.2998385 0.7473298 0.7075320
## V9
## 0.7029236fa_model = fa(EFA, nfactors = 3, fm = "ml", rotate = "promax")
factor_scores = factor.scores(EFA, fa_model)
EFA = cbind(EFA, factor_scores$scores)
names(EFA)## [1] "V1" "V2" "V3" "V4" "V5" "V6" "V7" "V8" "V9" "ML1" "ML3" "ML2"colnames(EFA)[colnames(EFA) == "ML1"] = "Faktor1"
colnames(EFA)[colnames(EFA) == "ML2"] = "Faktor2"
colnames(EFA)[colnames(EFA) == "ML3"] = "Faktor3"
names(EFA)## [1] "V1" "V2" "V3" "V4" "V5" "V6" "V7"
## [8] "V8" "V9" "Faktor1" "Faktor3" "Faktor2"Konfirmatorische Faktorenanalyse (CFA)
Hypothese: Die Zusammenhänge zwischen den Subtests Arithmetische Kompetenz, Figural-Induktives Denken, Wortbedeutung und Langzeitgedächtnis können mit einem Generalfaktor der Intelligenz erklärt werden.
# install.packages("lavaan")
# library(lavaan)
# Schritt 1: Modell definieren
CFA_model = "Generalfaktor =~ V1 + V2 + V3 + V4"
# Schritt 2: CFA durchführen
fit = cfa(CFA_model, data = CFA)
# Schritt 3: Ergebnisse anzeigen
summary(fit)## lavaan 0.6-18 ended normally after 28 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 8
##
## Number of observations 196
##
## Model Test User Model:
##
## Test statistic 2.356
## Degrees of freedom 2
## P-value (Chi-square) 0.308
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## Generalfaktor =~
## V1 1.000
## V2 1.231 0.320 3.848 0.000
## V3 1.110 0.293 3.790 0.000
## V4 0.965 0.262 3.679 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .V1 1.029 0.125 8.220 0.000
## .V2 0.644 0.111 5.776 0.000
## .V3 0.830 0.114 7.262 0.000
## .V4 0.787 0.100 7.833 0.000
## Generalfaktor 0.259 0.105 2.472 0.013# Erster Blick: Model Test User Model: Ergebnis n.s., d.h.: Die empirische Kovarianzmatrix entspricht in der Population der implizierten Kovariationsmatrix ("Model passt!")
summary(fit, fit.measures = TRUE)## lavaan 0.6-18 ended normally after 28 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 8
##
## Number of observations 196
##
## Model Test User Model:
##
## Test statistic 2.356
## Degrees of freedom 2
## P-value (Chi-square) 0.308
##
## Model Test Baseline Model:
##
## Test statistic 72.216
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.995
## Tucker-Lewis Index (TLI) 0.984
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1122.279
## Loglikelihood unrestricted model (H1) -1121.101
##
## Akaike (AIC) 2260.558
## Bayesian (BIC) 2286.783
## Sample-size adjusted Bayesian (SABIC) 2261.440
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.030
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.148
## P-value H_0: RMSEA <= 0.050 0.468
## P-value H_0: RMSEA >= 0.080 0.345
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.025
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## Generalfaktor =~
## V1 1.000
## V2 1.231 0.320 3.848 0.000
## V3 1.110 0.293 3.790 0.000
## V4 0.965 0.262 3.679 0.000
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .V1 1.029 0.125 8.220 0.000
## .V2 0.644 0.111 5.776 0.000
## .V3 0.830 0.114 7.262 0.000
## .V4 0.787 0.100 7.833 0.000
## Generalfaktor 0.259 0.105 2.472 0.013# Indizes CFI, RMSEA und SRMR sind relevant
## CFI = Comparative Fit Index
## RMSEA = Root Mean Square Error of Approximation
## SRMR = Standardized Root Mean Square Residual
summary(fit, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-18 ended normally after 28 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 8
##
## Number of observations 196
##
## Model Test User Model:
##
## Test statistic 2.356
## Degrees of freedom 2
## P-value (Chi-square) 0.308
##
## Model Test Baseline Model:
##
## Test statistic 72.216
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.995
## Tucker-Lewis Index (TLI) 0.984
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1122.279
## Loglikelihood unrestricted model (H1) -1121.101
##
## Akaike (AIC) 2260.558
## Bayesian (BIC) 2286.783
## Sample-size adjusted Bayesian (SABIC) 2261.440
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.030
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.148
## P-value H_0: RMSEA <= 0.050 0.468
## P-value H_0: RMSEA >= 0.080 0.345
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.025
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Generalfaktor =~
## V1 1.000 0.509 0.449
## V2 1.231 0.320 3.848 0.000 0.627 0.616
## V3 1.110 0.293 3.790 0.000 0.565 0.527
## V4 0.965 0.262 3.679 0.000 0.491 0.484
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .V1 1.029 0.125 8.220 0.000 1.029 0.799
## .V2 0.644 0.111 5.776 0.000 0.644 0.621
## .V3 0.830 0.114 7.262 0.000 0.830 0.722
## .V4 0.787 0.100 7.833 0.000 0.787 0.765
## Generalfaktor 0.259 0.105 2.472 0.013 1.000 1.000# install.packages("semPlot")
# library(semPlot)
semPaths(fit)
semPaths(fit, "std", edge.label.cex = 1, layout = "tree", whatLabels = "std")
names(CFA)## [1] "sex" "age" "V1" "V2" "V3" "V4"# str(CFA$sex)
CFA$sex = as.factor(CFA$sex)
# Drei Stufen der Messinvarianz: Konfigural, Metrisch, Skalare
# Konfigurale Messinvarianz
## "Passt das Modell in beiden Gruppen überhaupt, also sind dieselben Items überhaupt in denselben Faktoren enthalten?"
fit1 = cfa(CFA_model, data = CFA, group = "sex")
# Metrische Messinvarianz
## "Messen die Items den Faktor in beiden Gruppen gleich, also ist die Stärke mit der die Items auf einen Faktor laden, gleich?"
fit2 = cfa(CFA_model, data = CFA, group = "sex", group.equal = "loadings")
# Skalare Messinvarianz
## "Haben die Gruppen vergleichbare "Startpunkte" in der Messung?"
fit3 = cfa(CFA_model, data = CFA, group = "sex", group.equal = c("intercepts", "loadings"))
# Messinvarianzprüfung (Modellvergleich)
lavTestLRT(fit1, fit2, fit3)##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit1 4 2271.6 2350.3 5.0780
## fit2 7 2267.6 2336.4 7.0245 1.9465 0.00000 3 0.58359
## fit3 10 2271.4 2330.4 16.8817 9.8571 0.15272 3 0.01982 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Clusteranalyse
Es soll untersucht werden, ob sich Gruppen bzw. Typen an Personen identifizieren lassen, die sich hinsichtlich ihrer Giant-Three Persönlichkeitsprofile voneinander unterscheiden.
Big Three Model (Eysenck, 1994) Drei grundlegende Persönlichkeitseigenschaften: Extraversion, Neurotizismus und Psychotizismus (siehe Datensatz zur EFA)
CA$P = scale(CA$P)
CA$E = scale(CA$E)
CA$N = scale(CA$N)
head(CA)## # A tibble: 6 × 3
## P[,1] E[,1] N[,1]
## <dbl> <dbl> <dbl>
## 1 0.101 -0.260 -0.160
## 2 -1.50 -0.203 0.511
## 3 0.272 1.45 -1.30
## 4 -1.67 -0.317 -0.318
## 5 -2.07 -0.317 -0.633
## 6 0.272 -0.545 1.22Agglomeratives Hierarchisches Clustering
Berechnung der (Un-)Ähnlichkeit - Distanzmatrix
# Berechnung der paarweisen Distanzen zwischen den Zeilen
d = dist(CA, method = "euclidean")
as.matrix(d)[1:4, 1:4]## 1 2 3 4
## 1 0.000000 1.7360575 2.065952 1.7797259
## 2 1.736057 0.0000000 3.028084 0.8534846
## 3 2.065952 3.0280835 0.000000 2.8067161
## 4 1.779726 0.8534846 2.806716 0.0000000Durchführung des Clusteralgorithmus
Inkl. Linkage, typischerweise: Complete-Linkage-Methode (Maximum-Clustering) oder Ward-Methode
hc = hclust(d, method = "ward.D2")
# hc = hclust(d, method = "complete")
hc##
## Call:
## hclust(d = d, method = "ward.D2")
##
## Cluster method : ward.D2
## Distance : euclidean
## Number of objects: 299# install.packages("factoextra")
# library(factoextra)
fviz_dend(hc, cex = 0.2) + coord_flip()
Modellpassung
Validierung: Wie gut spiegelt das Dendrogramm die tatsächlichen Daten wider?
\(\rightarrow\) Korrelation zwischen cophenetischer Distanz und den ursprünglichen Distanzen
coph = cophenetic(hc) # Cophenetische Distanzen des Dendrogramms werden berechnet
cor(d, coph)## [1] 0.5399856# hc = hclust(d, method = "complete") # Complete-Linkage (Maximum Distance)
# coph = cophenetic(hc)
#
# cor(coph, d)
#
#
#
# hc = hclust(d, method = "single") # Single-Linkage (Minimum Distance)
# coph = cophenetic(hc)
#
# cor (coph, d)Nun: Zuweisung der Fälle zu Clustern (Tatsächliche Benennung der Cluster) \(\rightarrow\) Anzahl der Cluster k
Wahl der optimalen Anzahl von Clustern
Elbow-Methode
Berechnung der Gesamt-Within-Cluster Summe der Quadrate (WSS) (misst, wie ähnlich die Datenpunkte innerhalb eines Clusters sind) für verschiedene k-Werte (k = Anzahl der Cluster)
Identifikation des “Knick” (Elbow), bei dem die WSS abflacht
elbow_plot = fviz_nbclust(x = CA, FUN = hcut, method = c("wss"))
elbow_plot + geom_vline(xintercept = 3, linetype = 2, color = "black") 
# Manuelle Festlegung der Clusteranzahl durch "xintercept"
#
k_optimal = which.max(diff(diff(elbow_plot$data$y))) + 1
elbow_plot + geom_vline(xintercept = k_optimal, linetype = 2, color = "black")
Alternative: Silhouetten-Analyse
Misst die Qualität eines Clusterings anhand der Kohäsion (d.h. des Zusammenhalts innerhalb eines Clusters) und der Separation (Trennung der Cluster).
Silhouetten-Koeffizient von -1 bis +1: Werte nahe +1 deuten auf gut getrennte Cluster hin
fviz_nbclust(x = CA, FUN = hcut, method = c("silhouette"))
\(\rightarrow\) k = 2 ist der ideale Wert!
Zuweisung der Fälle zu den Clustern
# Dendrogramm nach k einfärben
fviz_dend(hc, cex = 0.2, k = 2, color_labels_by_k = TRUE) + coord_flip()
CA$cluster = cutree(hc, k = 2)
names(CA)## [1] "P" "E" "N" "cluster"Clusterbeschreibung und Interpretation
# Deskriptivstatistiken der Cluster
aggregate(CA[, c("P", "E", "N")], by = list(Cluster = CA$cluster), FUN = mean)## Cluster P E N
## 1 1 -0.06465722 0.3159969 -0.4086034
## 2 2 0.20017166 -0.9782917 1.2649914aggregate(CA[, c("P", "E", "N")], by = list(Cluster = CA$cluster), FUN = sd)## Cluster P E N
## 1 1 1.0333320 0.8833778 0.6945988
## 2 2 0.8650689 0.6453664 0.6961974Cluster 1:
Psychotizismus: Leicht unterdurchschnittlich, also tendenziell wenig impulsiv aggressiv.
Extraversion: Eher überdurchschnittlich, also eher gesellig, energisch und gesprächig
Neurotizismus: Eher unterdurchschnittlich, also emotional stabil, ruhig und weniger anfällig für negative Emotionen.
Mitglieder dieses Clusters sind tendenziell emotional stabil und gesellig, mit leicht unterdurchschnittlicher Impulsivität.
Cluster 2: Psychotizismus: Überdurchschnittlich, also tendenziell impulsiver/aggressiver.
Extraversion: Stark unterdurchschnittlich, also eher introvertiert, weniger gesellig und zurückhaltend.
Neurotizismus: Stark überdurchschnittlich, also emotional instabil, anfällig für Stress und negative Emotionen.
Mitglieder dieses Clusters sind eher introvertiert und emotional instabil, mit einer Tendenz zu impulsivem Verhalten.
ggplot(CA, aes(x = as.factor(cluster), y = P)) + geom_boxplot() + labs(title = "Clusterunterschiede: Psychotizismus")
ggplot(CA, aes(x = as.factor(cluster), y = E)) + geom_boxplot() + labs(title = "Clusterunterschiede: Extraversion")
ggplot(CA, aes(x = as.factor(cluster), y = N)) + geom_boxplot() + labs(title = "Clusterunterschiede: Neurotizismus")
t.test(P ~ cluster, data = CA)##
## Welch Two Sample t-test
##
## data: P by cluster
## t = -2.1641, df = 143.88, p-value = 0.03211
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -0.50671652 -0.02294124
## sample estimates:
## mean in group 1 mean in group 2
## -0.06465722 0.20017166# # ANOVA
# mod_P = lm(P ~ cluster, data = CA)
# anova(mod_P)
#
# # install.packages("emmeans")
# library(emmeans)
#
# CA$cluster = as.factor(CA$cluster)
#
# emmeans(mod_P, pairwise ~ cluster, adjust = "bonferroni")
t.test(E ~ cluster, data = CA)##
## Welch Two Sample t-test
##
## data: E by cluster
## t = 13.525, df = 166.06, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## 1.105344 1.483233
## sample estimates:
## mean in group 1 mean in group 2
## 0.3159969 -0.9782917t.test(N ~ cluster, data = CA)##
## Welch Two Sample t-test
##
## data: N by cluster
## t = -17.867, df = 121.72, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
## -1.859032 -1.488158
## sample estimates:
## mean in group 1 mean in group 2
## -0.4086034 1.2649914Strukturgleichungsmodelle
PoliticalDemocracy - Datensatz (Bollen, 1989) aus dem lavaan-Paket.
Latente Variablen (Faktoren), die im Modell spezifiziert werden sollen:
Demokratie im Jahr 1960 (dem60),
Demokratie im Jahr 1965 (dem65) und die
Industrielle Entwicklung im Jahr 1960 (ind60).
Diese werden durch manifeste Variablen (Indikatoren) gemessen:
dem60 wird durch die Indikatoren y1, y2, y3 und y4 repräsentiert,
dem65 wird durch die Indikatoren y5, y6, y7 und y8 repräsentiert,
ind60 wird durch die Indikatoren x1, x2 und x3 repräsentiert.
?PoliticalDemocracy
Hypothese: Die industrielle Entwicklung im Jahr 1960 (ind60) beeinflusst sowohl die Demokratie im Jahr 1960 (dem60) als auch die Demokratie im Jahr 1965 (dem65).
Es wird zudem angenommen, dass die Demokratie im Jahr 1960 (dem60) einen Einfluss auf die Demokratie im Jahr 1965 (dem65) hat.
Spezifizierung des Modells
Zwei Hauptkomponenten des SEMs:
Messmodell - Beschreibt die Beziehungen zwischen den latenten Variablen und ihren manifesten Indikatoren
Strukturmodell - Beschreibt die Beziehungen zwischen den latenten Variablen
Hinweis: In SEMs können die Residuen der Indikatoren miteinander korrelieren und müssen im Modell spezifiziert werden, um diesen zusätzlichen Zusammenhang im Modell zu berücksichtigen.
model = '
# Messmodell
ind60 =~ x1 + x2 + x3
dem60 =~ y1 + y2 + y3 + y4
dem65 =~ y5 + y6 + y7 + y8
# Strukturmodell
dem65 ~ ind60 + dem60
dem60 ~ ind60
# Residualkorrelationen
y1 ~~ y5
y2 ~~ y4 + y6
y3 ~~ y7
y4 ~~ y8
y6 ~~ y8
'Modell fitten
fit = sem(model, data = PoliticalDemocracy)
summary(fit, fit.measures = TRUE, standardized = TRUE, rsquare = TRUE)## lavaan 0.6-18 ended normally after 68 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 31
##
## Number of observations 75
##
## Model Test User Model:
##
## Test statistic 38.125
## Degrees of freedom 35
## P-value (Chi-square) 0.329
##
## Model Test Baseline Model:
##
## Test statistic 730.654
## Degrees of freedom 55
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.995
## Tucker-Lewis Index (TLI) 0.993
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1547.791
## Loglikelihood unrestricted model (H1) -1528.728
##
## Akaike (AIC) 3157.582
## Bayesian (BIC) 3229.424
## Sample-size adjusted Bayesian (SABIC) 3131.720
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.035
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.092
## P-value H_0: RMSEA <= 0.050 0.611
## P-value H_0: RMSEA >= 0.080 0.114
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.044
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ind60 =~
## x1 1.000 0.670 0.920
## x2 2.180 0.139 15.742 0.000 1.460 0.973
## x3 1.819 0.152 11.967 0.000 1.218 0.872
## dem60 =~
## y1 1.000 2.223 0.850
## y2 1.257 0.182 6.889 0.000 2.794 0.717
## y3 1.058 0.151 6.987 0.000 2.351 0.722
## y4 1.265 0.145 8.722 0.000 2.812 0.846
## dem65 =~
## y5 1.000 2.103 0.808
## y6 1.186 0.169 7.024 0.000 2.493 0.746
## y7 1.280 0.160 8.002 0.000 2.691 0.824
## y8 1.266 0.158 8.007 0.000 2.662 0.828
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## dem65 ~
## ind60 0.572 0.221 2.586 0.010 0.182 0.182
## dem60 0.837 0.098 8.514 0.000 0.885 0.885
## dem60 ~
## ind60 1.483 0.399 3.715 0.000 0.447 0.447
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .y1 ~~
## .y5 0.624 0.358 1.741 0.082 0.624 0.296
## .y2 ~~
## .y4 1.313 0.702 1.871 0.061 1.313 0.273
## .y6 2.153 0.734 2.934 0.003 2.153 0.356
## .y3 ~~
## .y7 0.795 0.608 1.308 0.191 0.795 0.191
## .y4 ~~
## .y8 0.348 0.442 0.787 0.431 0.348 0.109
## .y6 ~~
## .y8 1.356 0.568 2.386 0.017 1.356 0.338
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .x1 0.082 0.019 4.184 0.000 0.082 0.154
## .x2 0.120 0.070 1.718 0.086 0.120 0.053
## .x3 0.467 0.090 5.177 0.000 0.467 0.239
## .y1 1.891 0.444 4.256 0.000 1.891 0.277
## .y2 7.373 1.374 5.366 0.000 7.373 0.486
## .y3 5.067 0.952 5.324 0.000 5.067 0.478
## .y4 3.148 0.739 4.261 0.000 3.148 0.285
## .y5 2.351 0.480 4.895 0.000 2.351 0.347
## .y6 4.954 0.914 5.419 0.000 4.954 0.443
## .y7 3.431 0.713 4.814 0.000 3.431 0.322
## .y8 3.254 0.695 4.685 0.000 3.254 0.315
## ind60 0.448 0.087 5.173 0.000 1.000 1.000
## .dem60 3.956 0.921 4.295 0.000 0.800 0.800
## .dem65 0.172 0.215 0.803 0.422 0.039 0.039
##
## R-Square:
## Estimate
## x1 0.846
## x2 0.947
## x3 0.761
## y1 0.723
## y2 0.514
## y3 0.522
## y4 0.715
## y5 0.653
## y6 0.557
## y7 0.678
## y8 0.685
## dem60 0.200
## dem65 0.961semPaths(fit)
semPaths(fit,
edge.color = "black", # Pfade in schwarz
edge.label.cex = 0.8, # Größe der Labelschrift
edge.width = 0.5, # Einheitliche Liniendicke,
layout = "tree", # Layout (Alternativen: circle, spring)
sizeMan = 6, # Größe der manifesten Variablen
sizeLat = 8, # Größe der latenten Variablen
whatLabels = "std")
Anwendungsvoraussetzungen
pairs(PoliticalDemocracy[, 1:4]) # dem60
# pairs(subset(PoliticalDemocracy, select = c("y1", "y2", "y3", "y4")))
pairs(PoliticalDemocracy[, 5:8]) # dem65
pairs(PoliticalDemocracy[, 9:11]) # ind60
# Multikollinearität?
## Determinante der Korrelationsmatrix
cor_1 = cor(PoliticalDemocracy[, 1:4])
det(cor_1)## [1] 0.1197388cor_2 = cor(PoliticalDemocracy[, 5:8])
det(cor_2)## [1] 0.1016993cor_3 = cor(PoliticalDemocracy[, 9:11])
det(cor_3)## [1] 0.05381497# Stichprobengröße
nrow(PoliticalDemocracy)## [1] 75Mediation
Annahme: Es wird angenommen, dass der Zusammenhang zwischen Abhängigkeitskognitionen und der Ausprägung der Depressivität durch die negative Selbstbewertung vermittelt wird.
Hypothese im Sinne einer vollständigen Mediation
model <- ' # direct effect
BDI ~ c*ABK
# mediator
NSB ~ a*ABK
BDI ~ b*NSB
# indirect effect (a*b)
ab := a*b
# total effect
total := c + (a*b)
'
set.seed(123)
fit = sem(model, data = FIE, se = "bootstrap", bootstrap = 1000)
summary(fit, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-18 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 5
##
## Number of observations 191
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 176.345
## Degrees of freedom 3
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1353.544
## Loglikelihood unrestricted model (H1) -1353.544
##
## Akaike (AIC) 2717.088
## Bayesian (BIC) 2733.349
## Sample-size adjusted Bayesian (SABIC) 2717.511
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Bootstrap
## Number of requested bootstrap draws 1000
## Number of successful bootstrap draws 998
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## BDI ~
## ABK (c) 0.114 0.088 1.296 0.195 0.114 0.083
## NSB ~
## ABK (a) 0.507 0.082 6.161 0.000 0.507 0.420
## BDI ~
## NSB (b) 0.771 0.078 9.866 0.000 0.771 0.681
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .BDI 60.722 6.353 9.558 0.000 60.722 0.482
## .NSB 80.736 7.806 10.343 0.000 80.736 0.823
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ab 0.391 0.074 5.288 0.000 0.391 0.286
## total 0.505 0.091 5.571 0.000 0.505 0.369Der Zusammenhang zwischen Abhängigkeitskognitionen und der Ausprägung der Depressivität wird vollständig durch die negative Selbstbewertung vermittelt. Der direkte Effekte von Abhängigkeitskognitionen auf Depressivität ist nicht signifikant, während sowohl der Effekt von Abhängigkeitskognitionen auf die negative Selbstbewertung als auch der Effekt der negativen Selbstbewertung auf die Depressivität signifikant sind.
Praktische Implikationen: z.B. Therapeutische Ansätze: Interventionen sollten gezielt die negativen Selbstbewertungen ansprechen und verändern, da sie den entscheidenden Faktor darstellen.
Prävention: Menschen mit Abhängigkeitskognitionen könnten durch präventive Maßnahmen unterstützt werden, um negative Selbstbewertungen zu vermeiden und so das Risiko für Depressionen zu senken.